In [96]:
from __future__ import division, print_function

import os

# Third-party
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

# Custom
import gary.dynamics as gd
import gary.integrate as gi
import gary.io as io
import gary.potential as gp
from gary.units import galactic

from streammorphology import potential_registry, read_allfreqs

In [97]:
potential = potential_registry['via-lactea-log']

Read in initial conditions


In [40]:
w0 = np.loadtxt("/Users/adrian/projects/morphology/misc/progenitor_ICs.txt")
w0[:,3:6] = (w0[:,3:6]*u.km/u.s).decompose(galactic).value
classification = w0[:,6]
w0 = w0[:,:6]

Make a directory for freqmapping and save initial conditions with correct units


In [5]:
path = os.path.join("/Users/adrian/projects/morphology/output/ViaLacteaStreams/")
if not os.path.exists(path):
    os.mkdir(path)

In [8]:
w0_filename = os.path.join(path, "w0.npy")
np.save(w0_filename, w0)

In [41]:
t,all_orbits = potential.integrate_orbit(w0, dt=2., nsteps=25000, 
                                         Integrator=gi.DOPRI853Integrator)

In [42]:
per = np.sqrt(np.sum(all_orbits[:,:,:3]**2,axis=-1)).min(axis=0)
apo = np.sqrt(np.sum(all_orbits[:,:,:3]**2,axis=-1)).max(axis=0)


In [92]:
t1,w1 = potential.integrate_orbit(w0, dt=0.5, nsteps=100000, 
                                  Integrator=gi.DOPRI853Integrator)

In [93]:
naff = gd.NAFF(t1[:len(t1)//2+1])

In [94]:
i = 0

fs = [(w1[:len(t1)//2+1,i,j] + 1j*w1[:len(t1)//2+1,i,j+3]) for j in range(3)]
fxyz1,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)

fs = [(w1[len(t1)//2:,i,j] + 1j*w1[len(t1)//2:,i,j+3]) for j in range(3)]
fxyz2,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)

In [95]:
(fxyz2 - fxyz1) / fxyz1n


Out[95]:
array([  9.88617522e-06,  -1.00410444e-05,  -1.51150185e-06])


In [98]:
d = read_allfreqs(os.path.join(path, 'allfreqs.dat'), len(w0))
d.dtype.names


Out[98]:
('fxyz', 'fRphiz', 'dEmax', 'success', 'loop', 'dt', 'nsteps')

In [100]:
loop = d[d['loop']]
boxy = d[~d['loop']]

In [111]:
loop_nan = np.any(np.any(np.isnan(loop['fRphiz']), axis=-1), axis=-1)
boxy_nan = np.any(np.any(np.isnan(boxy['fxyz']), axis=-1), axis=-1)

In [116]:
loop[loop_nan]['dt'], loop[loop_nan]['nsteps']


Out[116]:
(array([ 0.69,  1.16,  1.86,  2.02,  3.56]),
 array([100000, 100000, 100000, 100000, 100000]))

In [117]:
boxy[boxy_nan]['dt'], boxy[boxy_nan]['nsteps']


Out[117]:
(array([ 1.33,  1.94,  2.58,  2.13,  2.66,  3.32,  2.74]),
 array([100000, 100000, 100000, 100000, 100000, 100000, 100000]))

Hm, ok, from the above, it looks like these orbits weren't integrated for long enough?


Plot the frequency diffusion rate vs. Emily's by-eye classification, colored by apocenter distance


In [118]:
boxy_fdiff = (boxy['fxyz'][:,1] - boxy['fxyz'][:,0]) / boxy['fxyz'][:,0] / (boxy['dt']*boxy['nsteps']/2.)[:,np.newaxis]
boxy_fdiff = np.max(boxy_fdiff, axis=1)
# boxy_fdiff = boxy_fdiff[np.isfinite(boxy_fdiff)]

plt.figure(figsize=(10,8))
plt.scatter(classification[~d['loop']], boxy_fdiff, c=apo[~d['loop']])
plt.colorbar()
plt.yscale('log')
plt.xlim(-1.,3.)
plt.ylim(1E-14, 1E-3)


Out[118]:
(1e-14, 0.001)

Ok, clearly I should only take apo < 50 kpc or something


In [145]:
ix = apo[~d['loop']] < 45

boxy_fdiff = (boxy['fxyz'][ix,1] - boxy['fxyz'][ix,0]) / boxy['fxyz'][ix,0] / (boxy['dt']*boxy['nsteps']/2.)[ix,np.newaxis]
boxy_fdiff = np.max(boxy_fdiff, axis=1)

loop_fdiff = (loop['fRphiz'][ix,1] - loop['fRphiz'][ix,0]) / loop['fRphiz'][ix,0] / (loop['dt']*loop['nsteps']/2.)[ix,np.newaxis]
loop_fdiff = np.max(loop_fdiff, axis=1)

# for coloring points
boxy_maxfreq = np.max(np.abs(boxy['fxyz'][ix,0]),axis=-1)
loop_maxfreq = np.max(np.abs(loop['fRphiz'][ix,0]),axis=-1)

fig,axes = plt.subplots(1,2,figsize=(12,6),sharex=True,sharey=True)

ax = axes[0]
c = ax.scatter(classification[~d['loop']][ix], boxy_fdiff, 
               marker='s', alpha=0.5, s=35, c='k')# c=per[~d['loop']][ix]) # c=boxy_maxfreq)
# fig.colorbar(c)
ax.set_yscale('log')
ax.set_xlim(-0.5,2.5)
ax.set_ylim(1E-14, 1E-3)
ax.xaxis.set_ticks([0,1,2])
ax.xaxis.set_ticklabels(['Normal','Fanned','Crazy'])
ax.set_title("Boxy", fontsize=20)

ax = axes[1]
c = ax.scatter(classification[d['loop']][ix], loop_fdiff,
               marker='s', alpha=0.5, s=35, c='k') # c=per[d['loop']][ix]) # c=loop_maxfreq)
# fig.colorbar(c)
ax.set_yscale('log')
ax.set_xlim(-0.5,2.5)
ax.set_ylim(1E-14, 1E-3)
ax.xaxis.set_ticks([0,1,2])
ax.xaxis.set_ticklabels(['Normal','Fanned','Crazy'])
ax.set_title("Loop", fontsize=20)


Out[145]:
<matplotlib.text.Text at 0x12bb330d0>

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [48]:
loop_fdiff = (loop['fRphiz'][:,1] - loop['fRphiz'][:,0]) / loop['fRphiz'][:,0] / (loop['dt']*loop['nsteps']/2.)[:,np.newaxis]
loop_fdiff = np.max(loop_fdiff, axis=1)
# loop_fdiff = loop_fdiff[np.isfinite(loop_fdiff)]

plt.figure(figsize=(10,8))
plt.scatter(classification[d['loop']], loop_fdiff, c=apo[d['loop']])
plt.colorbar()
plt.yscale('log')
plt.xlim(-1.,3.)
plt.ylim(1E-14, 1E-3)

plt.figure(figsize=(10,8))
plt.scatter(classification[d['loop']], loop_fdiff, c=per[d['loop']])
plt.colorbar()
plt.yscale('log')
plt.xlim(-1.,3.)
plt.ylim(1E-14, 1E-3)


Out[48]:
(1e-14, 0.001)

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [66]:
these_w0 = w0[w0[:,6] == 1.,:6][:10]

In [67]:
nsteps = 120000
t,w = potential.integrate_orbit(these_w0, dt=1., nsteps=nsteps, 
                                Integrator=gi.DOPRI853Integrator)

In [68]:
for i in range(w.shape[1]):
    fig = gd.plot_orbits(w[:nsteps//2], ix=i, marker='.', linestyle='none', alpha=0.1)
    fig = gd.plot_orbits(w[nsteps//2:], ix=i, marker='.', linestyle='none', alpha=0.1, axes=fig.axes)



In [69]:
for i in range(w.shape[1]):
    E = potential.total_energy(w[:,i,:3], w[:,i,3:6])
    dE = np.abs(E[1:] - E[0])
    print(dE.max())


3.60822483003e-16
4.71844785466e-16
5.27355936697e-16
7.91033905045e-16
1.12410081243e-15
4.16333634234e-16
1.27675647832e-15
4.99600361081e-16

In [62]:
naff = gd.NAFF(t[:len(t)//2+1])
fdiff = []
for i in range(w.shape[1]):
    fs = [(w[:len(t)//2+1,i,j] + 1j*w[:len(t)//2+1,i,j+3]) for j in range(3)]
    fxyz1,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)
    
    if (t.max() / np.abs(2*np.pi/fxyz1).max()) < 100:
        print("Skipping {}".format(i))
        continue
    
    fs = [(w[len(t)//2:,i,j] + 1j*w[len(t)//2:,i,j+3]) for j in range(3)]
    fxyz2,d,ixes = naff.find_fundamental_frequencies(fs, nintvec=15)
    
    # diffusion per orbit
    fdiff.append(np.abs((fxyz2 - fxyz1) / fxyz1) / (t.max() / np.abs(2*np.pi/fxyz1).max() / 2.))


Skipping 7

In [65]:
np.array(map(list, fdiff)).max(axis=1), np.array(map(list, fdiff)).max(axis=1).mean()


Out[65]:
(array([  9.65245684e-07,   6.90907826e-09,   1.07448038e-05,
          4.48579071e-07,   3.07767076e-06,   1.60719772e-06,
          4.32427600e-06]), 3.0249545940019157e-06)

Fractional change of 1E-3 per orbit for 'crazy'


In [63]:
t.max() / np.abs(2*np.pi/fxyz1).max() / 2.


Out[63]:
47.624549556280733

In [ ]: